bite <- function(T){
bite <- (1.67*10^-4) * T * (T- 2.3) * (32.0 - T)^(1/2)
bite <- ifelse(is.nan(bite), 0, bite)
return(bite)
}
transmit <- function(T){
-(2.94*10^-3) * T * (T - 11.3) * (T - 41.9)
}
landcover_avg <- readRDS(here::here("risk_maps",
"data",
"processed_data",
"landcover_avgs.RData"))
# Air temperature
t.test(filter(landcover_avg, landcover == "Agriculture")$air_temperature,
filter(landcover_avg, landcover == "Urban")$air_temperature,
paired = TRUE)
##
## Paired t-test
##
## data: filter(landcover_avg, landcover == "Agriculture")$air_temperature and filter(landcover_avg, landcover == "Urban")$air_temperature
## t = -42.019, df = 64, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -2.084814 -1.895572
## sample estimates:
## mean of the differences
## -1.990193
# biting rate
t.test(filter(landcover_avg, landcover == "Agriculture")$biting_rate,
filter(landcover_avg, landcover == "Urban")$biting_rate,
paired = TRUE)
##
## Paired t-test
##
## data: filter(landcover_avg, landcover == "Agriculture")$biting_rate and filter(landcover_avg, landcover == "Urban")$biting_rate
## t = 3.2143, df = 64, p-value = 0.00205
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.007554222 0.032364029
## sample estimates:
## mean of the differences
## 0.01995913
# transmission prob
t.test(filter(landcover_avg, landcover == "Agriculture")$transmission_prob,
filter(landcover_avg, landcover == "Urban")$transmission_prob,
paired = TRUE)
##
## Paired t-test
##
## data: filter(landcover_avg, landcover == "Agriculture")$transmission_prob and filter(landcover_avg, landcover == "Urban")$transmission_prob
## t = -5.0609, df = 64, p-value = 3.767e-06
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -1.5089083 -0.6548028
## sample estimates:
## mean of the differences
## -1.081856
all_pixels_integrate <- readRDS(here::here("risk_maps",
"data",
"processed_data",
"all_pixels_integrate.RData"))
coordinates(all_pixels_integrate) <- ~ x + y
gridded(all_pixels_integrate) <- TRUE
air_temp_raster <- raster(all_pixels_integrate, "air_temperature")
biting_raster <- raster(all_pixels_integrate, "biting_rate")
transmission_raster <- raster(all_pixels_integrate, "transmission_prob")
pal <- colorRampPalette(c("darkblue", "blue", "steelblue1", "white", "yellow", "orange", "red"))
plot(air_temp_raster,
col = pal(50),
# xlim=c(0,55),
# ylim=c(-10,50),
# zlim=limits,
xlab="Longitude",
ylab = "Latitude")
title(TeX("Air temperature ($^o$C)"))
plot(biting_raster,
col = pal(50),
# xlim=c(0,55),
# ylim=c(-10,50),
# zlim=limits,
xlab="Longitude",
ylab = "Latitude")
title(TeX("Biting rate ($day^{-1}$)"))
plot(transmission_raster,
col = pal(50),
# xlim=c(0,55),
# ylim=c(-10,50),
# zlim=limits,
xlab="Longitude",
ylab = "Latitude")
title(TeX("Transmission probability (%)"))
Maps of errors integrated
get_error <- function(temp_raster, risk_function){
max_error_by_pixel <- function(temperature){
if (is.na(temperature)){
return(NA)
} else {
lower_bound <- temperature - (STANDARD_ERROR/(65-1))*1.96
upper_bound <- temperature + (STANDARD_ERROR/(65-1))*1.96
max <- optimize(risk_function, interval=c(lower_bound, upper_bound), maximum=TRUE)$objective
min <- optimize(risk_function, interval=c(lower_bound, upper_bound), maximum=FALSE)$objective
return(max - min)
}
}
values <- raster::values(temp_raster)
errors <- sapply(values, max_error_by_pixel)
values(temp_raster) <- errors
return(temp_raster)
}
plot(get_error(air_temp_raster, bite),
col = pal(50),
# xlim=c(0,55),
# ylim=c(-10,50),
# zlim=limits,
xlab="Longitude",
ylab = "Latitude")
title(TeX("Biting rate uncertainty ($day^{-1}$)"))
plot(get_error(air_temp_raster, transmit),
col = pal(50),
# xlim=c(0,55),
# ylim=c(-10,50),
# zlim=limits,
xlab="Longitude",
ylab = "Latitude")
title(TeX("Transmission probability uncertainty (%)"))
## NULL
Because each ECOSTRESS image has a slightly different grid, it is difficult to match the pixels to each other. We do this by resampling all images to the first image.
pixels <- readRDS(here::here("risk_maps",
"data",
"processed_data",
"all_pixels_integrate_w_landcover.RData")) %>%
filter(!is.na(air_temperature))
ag_types <- c("Vegetable", "Orchard", "Fruit", "Field Crop", "Uncultivated", "Urban")
pixels$landcover <- factor(pixels$landcover, levels = c("Urban", "Uncultivated", "Vegetable", "Fruit", "Orchard", "Field Crop", "Other"))
filter(pixels, landcover %in% ag_types) %>%
ggplot(aes(x = air_temperature, col= landcover, alpha = 0.2)) +
geom_density() +
guides(alpha = FALSE, col = guide_legend(title = "Crop category")) +
xlab(TeX("Air temperature ($^o$C)")) +
ylab("Density") +
xlim(22.5, 29.2) +
# labs(title = "Average air temperature distribution") +
theme_bw() +
theme(legend.position = c(0.15, 0.7), legend.title = element_blank())
filter(pixels, landcover %in% ag_types) %>%
ggplot(aes(x = transmission_prob, col= landcover, alpha = 0.2)) +
geom_density() +
guides(alpha = FALSE, col = guide_legend(title = "Crop category")) +
xlab("Transmission probability (%)") +
ylab("Density") +
xlim(13, 17.3) +
# labs(title = "Average transmission probability distribution") +
theme_bw() +
guides(color = FALSE)
filter(pixels, landcover %in% ag_types) %>%
ggplot(aes(x = biting_rate, col= landcover, alpha = 0.2)) +
geom_density() +
guides(alpha = FALSE, col = guide_legend(title = "Crop category")) +
xlab(TeX("Biting rate ($day^{-1}$)")) +
ylab("Density") +
xlim(.11, .225) +
# labs(title = "Average biting rate distribution") +
theme_bw() +
guides(color = FALSE)
# original submission: Instead of just looking at a single daytime frame's pixels, look at all measurements from 11 to 16h (15 frames)
# day_pixels <- readRDS(here::here("risk_maps",
# "data",
# "processed_data",
# "day_pixels.RData"))
#
# ag_types <- c("Vegetable", "Orchard", "Fruit", "Field Crop", "Uncultivated")
# day_pixels$landcover <- factor(day_pixels$landcover, levels = c("Urban", "Uncultivated", "Vegetable", "Fruit", "Orchard", "Field Crop", "Other"))
#
# filter(day_pixels, landcover %in% ag_types) %>%
# ggplot(aes(x = air_temperature, col= landcover, alpha = 0.2)) +
# geom_density() +
# guides(alpha = FALSE, col = guide_legend(title = "Crop category")) +
# xlab(TeX("Air temperature ($^o$C)")) +
# ylab("Density") +
# xlim(25, 36.7) +
# labs(title = "Daytime air temperature distribution") +
# theme_bw()
#
# filter(day_pixels, landcover %in% ag_types) %>%
# ggplot(aes(x = transmission_prob, col= landcover, alpha = 0.2)) +
# geom_density() +
# guides(alpha = FALSE, col = guide_legend(title = "Crop category")) +
# xlab("Transmission probability (%)") +
# ylab("Density") +
# xlim(14.5, 20) +
# labs(title = "Daytime transmission probability distribution") +
# theme_bw()
f_air_temp <- aov(air_temperature ~ landcover, data = filter(pixels, landcover %in% ag_types))
f_trans_prob <- aov(transmission_prob ~ landcover, data = filter(pixels, landcover %in% ag_types))
summary(f_air_temp)
## Df Sum Sq Mean Sq F value Pr(>F)
## landcover 5 333260 66652 120990 <2e-16 ***
## Residuals 448948 247321 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(f_trans_prob)
## Df Sum Sq Mean Sq F value Pr(>F)
## landcover 5 87220 17444 145652 <2e-16 ***
## Residuals 448948 53768 0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
####(space averages)
All <-filter(landcover_avg, landcover == "All")
bite <- function(T){
bite <- (1.67*10^-4) * T * (T- 2.3) * (32.0 - T)^(1/2)
bite[is.nan(bite)] <- 0
return(bite)
}
transmit <- function(T){
-(2.94*10^-3) * T * (T - 11.3) * (T - 41.9)
}
All$avg_bite <- All$air_temperature %>% bite
All$avg_transmit <- All$air_temperature %>% transmit
t.test(x = All$biting_rate, y = All$avg_bite, paired = TRUE)
##
## Paired t-test
##
## data: All$biting_rate and All$avg_bite
## t = 2.5515, df = 64, p-value = 0.01313
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.001369638 0.011251819
## sample estimates:
## mean of the differences
## 0.006310729
t.test(x = All$transmission_prob, y = All$avg_transmit, paired = TRUE)
##
## Paired t-test
##
## data: All$transmission_prob and All$avg_transmit
## t = -7.7387, df = 64, p-value = 9.314e-11
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.3090675 -0.1822379
## sample estimates:
## mean of the differences
## -0.2456527
####(time averages)
all_pixels_integrate <- readRDS(here::here("risk_maps",
"data",
"processed_data",
"all_pixels_integrate.RData")) %>% filter(!is.na(air_temperature)) #there are 3751 pixels on the edge that do not have values for every image, and these we remove.
all_pixels_integrate$avg_bite <- all_pixels_integrate$air_temperature %>% bite
all_pixels_integrate$avg_transmit <- all_pixels_integrate$air_temperature %>% transmit
t.test(x = all_pixels_integrate$biting_rate, y = all_pixels_integrate$avg_bite, paired = TRUE)
##
## Paired t-test
##
## data: all_pixels_integrate$biting_rate and all_pixels_integrate$avg_bite
## t = -4384.6, df = 964748, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.09894369 -0.09885528
## sample estimates:
## mean of the differences
## -0.09889949
t.test(x = all_pixels_integrate$transmission_prob, y = all_pixels_integrate$avg_transmit, paired = TRUE)
##
## Paired t-test
##
## data: all_pixels_integrate$transmission_prob and all_pixels_integrate$avg_transmit
## t = -3129.6, df = 964748, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -3.069540 -3.065697
## sample estimates:
## mean of the differences
## -3.067618
all_pixels_integrate <- readRDS(here::here("risk_maps",
"data",
"processed_data",
"all_pixels_integrate.RData")) %>% filter(!is.na(air_temperature)) #there are 3751 pixels on the edge that do not have values for every image, and these we remove.
avg_bite <- all_pixels_integrate$air_temperature %>% mean(na.rm = TRUE) %>% bite
avg_transmit <- all_pixels_integrate$air_temperature %>% mean(na.rm = TRUE) %>% transmit
t.test(x = all_pixels_integrate$biting_rate, mu = avg_bite)
##
## One Sample t-test
##
## data: all_pixels_integrate$biting_rate
## t = -4372.7, df = 964748, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0.2512659
## 95 percent confidence interval:
## 0.1492290 0.1493204
## sample estimates:
## mean of x
## 0.1492747
t.test(x = all_pixels_integrate$transmission_prob, mu = avg_transmit)
##
## One Sample t-test
##
## data: all_pixels_integrate$transmission_prob
## t = -5128.9, df = 964748, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 18.21251
## 95 percent confidence interval:
## 15.04580 15.04822
## sample estimates:
## mean of x
## 15.04701
bitemax <- optimize(bite, interval=c(0, 40), maximum=TRUE)$objective
transmax <- optimize(transmit, interval=c(0, 40), maximum=TRUE)$objective
ggplot(All) + geom_point(aes(x = avg_bite - biting_rate, y = hour(dt) + minute(dt)/60, color = avg_bite))
ggplot(All) + geom_point(aes(x = avg_transmit - transmission_prob, y = hour(dt) + minute(dt)/60, color = avg_transmit))
lstfiles <- list.files(path = here::here("risk_maps",
"data",
"raw_data",
"ECOSTRESS",
"all_ims",
"LST"), pattern = "tif", recursive = TRUE, full.names = FALSE)
getdate <- function(file){
nums <- str_extract(file, regex('(?<=doy)[0-9]+'))
year <- substr(nums, 1, 4)
doy <- substr(nums, 5, 7) %>% as.numeric()
hhmmss <- substr(nums, 8, 13)
date <- as.Date(doy - 1, origin = paste0(year, "-01-01"))
dt <- ymd_hms(paste(date, hhmmss), tz = "UTC") %>% with_tz("America/Los_Angeles")
return(dt)
}
all_dates <- do.call(c, lapply(lstfiles, getdate))
files <- list.files(path = here::here("risk_maps",
"data",
"raw_data",
"ECOSTRESS",
"filtered_ims"),
pattern = "tif",
recursive = TRUE,
full.names = FALSE)
getfiltdate <- function(file){
date <- ymd(substr(file, 14, 23))
hhmmss <- str_extract(file, regex('[0-9]{2}:{1}[0-9]{2}:{1}[0-9]{2}'))
dt <- ymd_hms(paste(date, hhmmss), tz = "America/Los_Angeles")
return(dt)
}
filtered_dates <- do.call(c, lapply(files, getfiltdate))
data.frame(dt = all_dates, included = ifelse(all_dates %in% filtered_dates, "Included", "Not included")) %>%
ggplot() +
geom_point(aes(x = month(dt) + day(dt)/30 + hour(dt)/(24*30), y = hour(dt) + minute(dt)/60, shape = as.factor(year(dt)), color = included)) +
scale_color_manual(values = c("red", "lightblue2")) +
xlab("Month") +
ylab("Hour") +
theme_bw() +
theme(legend.title = element_blank())